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Abstract 

We carry out 3-D numerical simulations of the dynamical instability in rapidly 
rotating stars initially modeled as polytropes with n = 1.5, 1.0, and 0.5. The 
calculations are done with a SPH code using Newtonian gravity, and the grav- 
itational radiation is calculated in the quadrupole limit. All models develop 
the global m = 2 bar mode, with mass and angular momentum being shed 
from the ends of the bar in two trailing spiral arms. The models then un- 
dergo successive episodes of core recontraction and spiral arm ejection, with 
the number of these episodes increasing as n decreases: this results in longer- 
lived gravitational wave signals for stiffer models. This instability may operate 
in a stellar core that has expended its nuclear fuel and is prevented from fur- 
ther collapse due to centrifugal forces. The actual values of the gravitational 
radiation amplitudes and frequencies depend sensitively on the radius of the 

star i?eq at which the instability develops. 

PACS numbers: 04.30.Db, 04.80.Nn, 95.55.Ym, 97.60.-s 



1 



Typeset using REVTJ5X 



I. INTRODUCTION 



The direct detection of gravitational radiation from astrophysical sources is one of the 
greatest challenges of our day. With interferometers such as LIGO [[I], VIRGO 0, and 
GEO [[| under construction, and a new generation of spherical resonant mass detectors 
under study [|]|| , the detailed modeling of these sources takes a high priority. 

One promising source of gravitational waves is the development of rotational instabilities 
in dense stellar cores or compact objects ||. Consider, for example, a rapidly rotating stellar 
core that has expended its nuclear fuel and is unable to collapse to neutron star size due to 
centrifugal forces. If such an object underwent a rotational instability, it could possibly shed 
enough angular momentum to allow collapse to a supernova ]7|||. Alternatively, a neutron 
star that is spun up by accretion of mass from a binary companion may reach fast enough 
rotation rates to go unstable [p|,|10| . 



Global rotational instabilities in fluids arise from nonradial "toroidal" modes e ±m< ^, 
where m = 2 is known as the "bar mode" . It is convenient to parametrize them by 

[3 = T rot /\W\, (1) 

where T rot is the rotational kinetic energy and W is the gravitational potential energy 



TT| , p~2| . p~3|] . We focus on the bar instability since it is expected to be the fastest growing 



mode. This instability can occur under two different physical mechanisms. The dynamical 
bar instability is driven by Newtonian hydrodynamics and gravity. It operates for fairly 
large values of the stability parameter (3 > (3d and develops on a timescale of approximately 
one bar rotation period. In contrast, the secular instability arises from dissipative processes 
such as gravitational radiation reaction and viscosity. It occurs for (3 S < (3 < (3d and de- 
velops on a timescale of several rotation periods or longer [JTIJ. For the constant density, 



incompressible, uniformly rotating Maclaurin spheroids we have (3 S ~ 0.14 and (3d ~ 0.27. 
In the case of differentially rotating fluids with a polytropic equation of state 

P = Kp T = Kp 1+1/n , (2) 



where n is the polytropic index and K is a constant that depends on the entropy, early 
studies indicated that the secular and dynamical bar instabilities should occur at about 
these same values of (3 [|l^,|l^,|14],|Uj . More recent work [TJJ shows that both the angular 



momentum distribution and, to a lesser degree, the polytropic index affect the value of f3 s 
at which the m = 2 secular instability sets in. For the dynamical bar instability Pickett, 



et al. |L7| demonstrate that, for n = 1.5 polytropes, the m = 2 dynamical stability limit 
fid ~ 0.27 is valid for centrally condensed initial angular momentum distributions that are 
similar to those of Maclaurin spheroids. However, for angular momentum distributions with 
somewhat extended disk-like regions, both one- and two-armed spiral instablities appear at 
considerably lower values of j3. 

The work presented in this paper is part of a research program aimed at calculating 
the gravitational radiation produced when a rapidly rotating stellar core undergoes the 
dynamical bar instability. These studies are carried out using 3-D numerical simulations. 
The gravitational field is purely Newtonian, and the gravitational radiation produced is 
calculated using the quadrupole approximation; the back reaction of the radiation on the 
system is not included. The rapidly rotating cores are initially modeled as polytropes with 
(3 ~ 0.3, which is just above the dynamical stability limit. In order to sustain such high 
rotational kinetic energy they must be rotating differentially fll|Jl8|| . Such objects could 
form, for example, when the cores of massive stars collapse on a dynamical time scale. 

Much of the previous work in this area has concentrated on polytropes with n = 1.5. 
This case has been investigated by Tohline and collaborators [O,p0,pl[] and more recently 



by Pickett, et al. [[L7[ in the context of star formation. These studies primarily use an 
Eulerian code that imposes the polytropic equation of state @ throughout the evolution 
instead of solving an energy equation; thus, the entropy generation by shocks during the 
later stages of evolution is not taken into account. The work of Houser, Centrella, and 



Smith was the first to model the fluid using an energy equation and to calculate the 
resulting gravitational radiation; these calculations were carried out using the Lagrangian 
smooth particle hydrodynamics (SPH) method. (New |23] has recently performed similar 



calculations using an improved version of Tohline's Eulerian code.) Smith, Houser, and 
Centrella [23[] also updated the earlier work of Ref. |2(J by carrying out a detailed comparison 
study of this model with two different 3-D codes, one using Eulerian techniques and the other 
based on SPH. 

In this paper we extend our calculations to objects having stiff polytropic equations of 
state using SPH |25[| . We use (3 ~ 0.3 and consider the cases n = 1.5, 1.0, and 0.5, which 



correspond to V = 5/3,2, and 3. Previously, Williams and Tohline [£l[] studied the initial 
development of the bar instability in similar models with n = 0.8, 1.0, 1.3, 1.5, and 1.8; longer 
evolutions were carried out in Ref. |26j for the cases n = 0.8 and n = 1.8. Their work was 
done with an Eulerian code in cylindrical coordinates (vj, <p, z, ) that used the diffusive donor 
cell advection method and imposed the polytropic equation of state (§) throughout the runs. 
In addition, they modeled only the region < <fi < tt in the angular coordinate, so that only 
even toroidal modes were represented. Our simulations do not suffer from these restrictions, 
and include the calculation of the gravitational radiation. 

This paper is organized as follows. Section [TJ contains a brief description of the techniques 
used in our simulations. The construction of initial axisymmetric models with f3 ~ 0.3 is 
discussed in Sec. |TT| and the dynamical evolution of the models is presented in Sec. [TV| . 
Analysis of the instabilities using Fourier components is given in Sec. [V| and the gravitational 



radiation produced by the models is presented in Sec. [Vl|. The paper concludes with a 
summary and discussion of our results in Sec. [VI]] . 

II. SIMULATION TECHNIQUES 

Detailed descriptions of the basic techniques used to carry out these simulations have 
been presented previously in Refs. and ||27|| . We therefore give a only brief description of 



these methods in this section, and refer the reader to the literature for further information. 

SPH is a Lagrangian method in which the fluid is modeled as a collection of fluid elements 
of finite extent described by a smoothing kernel We have used the implementation of 
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SPH by Hernquist & Katz |29| known as TREESPH, which allows variable smoothing lengths 
and individual particle timesteps. For the runs discussed in this paper, we smooth over 
N$ = 64 neighbors for kernel interpolation. Shocks are handled using an artificial viscosity 
modified by the curl of the velocity field, with the user-specified coefficients of the linear 
and quadratic terms taking the values a = 0.25 and /?av — 1-0; see Refs. [ P§| , F^| , pT| f° r m ore 



details. The gravitational forces in this code are purely Newtonian, and are calculated using 
a hierarchical tree method optimized for vector computers ||30|| . This leads to a significant 



gain in efficiency and allows the use of larger numbers of particles than would be possible 
with methods that simply sum over all possible pairs of particles. 

We calculate the gravitational radiation produced by the instabilities using the 



quadrupole approximation, which is valid for nearly Newtonian sources |3T]. The reduced 
(i.e., traceless) quadrupole moment of the source is given by 

hj = J P (xiXj ~ \5ijr 2 ) d 3 r, (3) 

where 2, j = 1, 2, 3 are spatial indices and r = (x 2 + y 2 + z 2 ) 1 ^ 2 is the distance to the source. 
For an observer situated on the axis at 9 = 0, = of a spherical coordinate system with its 
origin located at the center of mass of the source, the amplitude of the gravitational waves 
for the two polarization states takes the simple form 

G 1 

h + = - A -{hx-iyy)i ( 4 ) 
q 2 

= ^- f a>w, (5) 

where an overdot indicates a time derivative d/dt. The gravitational wave luminosity is 
given by 

and the rate at which angular momentum is lost through gravitational radiation is 

dJj _2G /r(2) f (3)\ (7 s 

dt ~ 5 c 5 Vim 1 km/ ■ v) 



Here there is an implied sum on repeated indices, the superscript (3) indicates the third 
time derivative, and the angle brackets indicate an average over several wave periods. Since 
such averaging is not well-defined for these burst sources, we display instead the unaveraged 
quantities (G/Sc^fgMf and (2G/5c 5 ) 

^ijkljm^km below. The energy emitted as gravita- 
tional radiation is 



AE = J Ldt (8) 
and the angular momentum carried away by the waves is 

AJi = J (dJi/dt)dt. (9) 
Finally, the gravitational wave energy spectrum dE/df, which gives the energy emitted as 



gravitational radiation per unit frequency interval, takes the form |32 



W = ^^^W + IM/)I 2 », (io) 



where h(f) is the Fourier transform of h{t). The double angle brackets in Eq. ( |I0D denote 
an average over all source angles. 

We calculate the reduced quadrupole moment fy and its derivatives using the methods 



developed in Ref. ||27|| . In particular, particle positions, velocities, and accelerations already 
present in the code are used to obtain 1^ and ly, yielding expressions similar to those in 
Ref. |33| . This results in waveforms that are very smooth functions of time and require no 
filtering or smoothing to remove numerical noise. However, the luminosity L and angular 
momentum lost by gravitational radiation dJi/dt do contain the time derivative of the par- 
ticle acceleration; this is taken numerically and therefore introduces some noise. To remove 
this noise, we smooth the luminosity data using simple averaging over a fixed time interval 
of O-lto centered on each point. Here, £d is the dynamical time for a spherical star of mass 
M and (equatorial) radius R eq and is defined by 

/ R 3 \ 1/2 



In general this procedure produces very smooth luminosity profiles and makes a negligi- 
ble change in the integrated luminosity AE, which gives the energy emitted as gravitational 
radiation. The profiles of dJi/dt are not smoothed. 



III. INITIAL MODELS 

The initial conditions for our simulations consist of rotating axisymmetric equilibrium 
fluid models having (3 ~ 0.3 and polytropic index n. The dynamical bar instability then 
grows from nonaxisymmetric perturbations due to particle discreteness in the models. We 
use a two-step procedure to generate the initial models. First, a self-consistent field (SCF) 
method is used to produce an equilibrium model on a grid. Then, a particle fit to the SCF 
model is performed to generate initial data for TREESPH. In this section, we describe the 
construction of these equilibrium models. 



A. Self-Consistent Field Method 



The first step is to use the SCF method ( see also p5| , |36| , p7| ) to generate axisym- 
metric equilibrium models. The SCF procedure derives from an integral formulation of the 
equations of hydrodynamic equilibrium which automatically incorporates the boundary con- 
ditions. We use cylindrical coordinates (w, z) and a uniformly- zoned grid of N m radial and 
N z axial zones. An initial "guess" density distribution p(w, z) is given, and the gravitational 
potential is calculated using a Legendre polynomial expansion to solve Poisson's equation 
|37| . A rotation law in which angular momentum is constant on cylinders is specified. This 
takes the general form j(m) = j(m(tz7)), where j(m) is the specific angular momentum and 



m{w) here denotes the dimensionless mass fraction interior to the cylinder of radius w p5 



Following the convention of earlier work (e.g. Refs. |L^,^(],^1],^,[!(J ) we use the rotation 
law for the uniformly rotating, constant density Maclaurin spheroids, 

s 5 J 



nm } = -- 



(l-m) 2/3 , (12) 



where J is the total angular momentum and M is the total mass. Applying this rotation law 
to polytropes, which do not have constant density, produces differentially rotating models. 
The rotation law ([X2|) is used to calculate a rotational potential, which is then used with the 
gravitational potential to compute an improved density distribution. This process is then 
repeated, iterating until convergence is achieved. 

Rotation causes the resulting models to be flattened, so that R p < R eq , where R p is 
the polar radius and R eci is the equatorial radius. The freely specifiable quantities in this 
method are the dimensionless form of the rotation law h(m) = (M/ J)j(m), n, i? eq , and 
the axis ratio R p /R eq . Upon convergence to a solution of the equations of hydrodynamic 
equilibrium, this procedure gives the density p(va,z), the angular velocity fl(w), M, J, and 
(3. To get a dimensional model, we specify the entropy, which is given by the constant K in 
the polytropic equation of state (0), and the maximum density. 

One measure of the accuracy of the initial equilibrium models comes from the virial 



theorem. For a fluid system this gives [ 1 1 



2T + W + 311 = 0, (13) 

where T is the total kinetic energy and II = / PdV is the volume integral of the pressure. 
Using this we define the virial relation VR by |37 



VR 



2T + W + 311 



(14) 



W 

Using the SCF method, we generated initial models for the cases n = 1.5, 1.0, and 0.5; 
we refer to these as SCF1, SCF2, and SCF3, respectively. The parameters of these initial 
models are given in Table |I[ Notice that a finer grid was used for the n = 1.5 model SCF1; 
this was done because it is more centrally condensed than the other two cases. 



B. Generation of Particle Models 



Once the SCF equilibrium model has been produced on a cylindrical grid, it must be 
transformed into a form readable by TREESPH. This is done by performing a particle 
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fit to the density profile p(w, z) given by the SCF method. The simplest technique for 
doing this randomly distributes particles within the probability distribution p(za, z) using 
a "rejection" method p8| , p?| ; this technique was used to generate the initial conditions for 



the runs performed in Ref. |[22|| . Due to the underlying Poisson distribution of particles, 
these initial models contain relatively large positive and negative density fluctuations. The 
resulting internal noise causes spurious entropy generation and partially masks the signal 
of the bar mode instability. Because of these problems, we developed methods to produce 
"colder" initial particle models with less internal noise |24|,f2g]. 

To generate these cold models, we use equipotential surfaces to determine the function 
M($), which is the mass interior to an equipotential. This information then allows us to 
obtain the desired physical properties of the SCF model by reordering the particles from 
a physically simpler model. This simpler model is created by placing particles within the 
known stellar boundary to obtain a uniform density particle distribution. Then, using a 
chosen set of equipotential surfaces, the mass interior to these surfaces M($) is computed in 
both the SCF and uniform density models. The actual number of surfaces used is taken to 
equal the number of zones in the w direction used to generate the SCF model, N pot = N m ; 
c.f. Table |. By a direct comparison between the resulting SCF and uniform density mass 
functions, a systematic contraction, or repositioning, of particle positions from their original 
locations in the uniform density model can be performed. This results in a particle model 
which realistically reproduces the SCF density profile, and does not suffer from the density 
fluctuations found in the random particle method . 



To implement this procedure, we use the rotational ($ ro t) and gravitational ($ gr(M) ) 
potentials, which are natural by-products of the SCF method, to define the total surface 
equipotential $ SOT /- A uniform 3-D Cartesian grid centered on the star is then created inside 
a cube having length 2R eq . Particles are placed at each of the grid nodes, and a particle is 
accepted into the model if it lies inside the boundary of the star, producing a 3-D uniform 
density particle representation with the exact physical shape of the SCF stellar model. This 
uniform distribution must now be transformed into the more centrally-condensed polytropic 
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model given by the SCF method. In practice we do this by systematically "contracting" the 
particle positions in the original uniform density model along radial vectors (i.e., moving 
them toward the center). Comparison of the mass functions M($) for the two models 
tells how to move the particles along radial rays to achieve the SCF density distribution. 
Since the SCF model is more centrally condensed than the uniform density distribution, 
the particles are moved toward the stellar center to their new positions. This repositioning 
of particles then reproduces the SCF density distribution. When the contraction process 
is completed, each particle is assigned an angular velocity Q to reproduce the rotation law 
given in Eq. ([12]). 

Using this contraction method, we generated initial models with (3 ~ 0.3 and n = 1.5 
using particle numbers in the range N ~ 2000 — 32, 000. These SPH models form the basis 
of the comparison study with Eulerian methods reported in Ref. p4| |. Overall, these models 



do not suffer from the large density fluctuations present in the randomly generated models 
and, when evolved using TREESPH, better reproduce the basic features of the bar mode 
instability. 

To generate the initial models for the runs presented here, we incorporated several im- 
provements to this method. First of all, we added an iterative procedure to the contraction 
process. The initial repositioning of particles is identical to that presented above. However, 
once the uniform density model has been contracted, the mass interior to the equipotential 
surfaces is recalculated for the particle model. The radial contraction is again applied to 
all particles, and the process is iterated until the difference between the initial and final 
positions of all particles is less than a given tolerance, here chosen to be typically ~ O.SiV" 1 . 
We also modified the initial placement of the particles within the uniform density model. A 
box with extent ±R eq in the x and y directions and ±R P in the z direction was used instead 
of a cubical box. Squeezing the particles along the z direction to fit within the range of R p 
should initially position them closer to their final equilibrium positions, thus making the 
contraction method more efficient. Also, to eliminate systematic errors due to contraction, 
the particle planes of constant z were displaced in the w direction by ±1/4 of the inter- 
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particle spacing for even and odd z-planes, respectively. Overall, this iterated contraction 
method produces particle models that better reproduce the SCF density distribution. 

Using this iterated contraction procedure, we constructed three models with (3 w 0.3 
using N ~ 16000 for different values of the polytropic index n. We refer to these as Run 1 
(n = 1.5), Run 2 (n = 1.0), and Run 3 (n = 0.5); see Table [TT| Fig. [I] shows the normalized 
equatorial plane density for these initial models. Notice that the density values calculated 
at the particle positions (shown as dots in the figure) match the SCF profiles (shown as 
solid lines) with very little scatter. The angular velocity profiles for all particles are shown 
in Fig. § and also reproduce the SCF values with very little scatter. 

Rigorously, this repositioning of particles should be carried out by following normal, 
rather than radial, vectors. We have used radial contraction here for simplicity and compu- 
tational speed. For spherical systems, the normal and radial vectors coincide. However, as 
the equilibrium model becomes increasingly oblate due to rotation, contracting the particles 
along their radial vectors becomes less accurate. In practice, we find the radial contraction 
method models the equatorial plane density well, as shown in Fig. [I]. However, under- 
densities are observed in the regions around the rotation axis for increasing values of \z\. 
When the models are evolved forward in time, this causes a slight redistribution of mass 
and angular momentum in the inner regions. Based on comparisons of the n = 1.5 singly 



contracted models with Eulerian runs in Ref. [24}| , we do not believe that this small adjust- 
ment of the initial models significantly affects the evolution of the bar mode instability. The 
improved iterated contraction method reduces this under-density somewhat for the n = 1.5 
case. As n decreases, the polytrope becomes less centrally condensed and hence is closer to 
the uniform density model prior to contraction. As a result the under-densities decrease for 
Run 2, and almost disappear for Run 3. 
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IV. DYNAMICAL EVOLUTION 



The initial particle models for Runs 1, 2, and 3 generated using the iterated contraction 
method were evolved in time using TREESPH. The case n — 1.5 has already been the 
subject of our detailed comparison study using Eulerian and SPH methods reported in Ref. 
||24|| . Run 1 comprises the evolution of an improved initial model with n = 1.5 and shows 
the same general behavior seen in that previous study; it is included here as a benchmark 
for comparison with the stiffer equations of state. 

The dynamical evolution of these models is displayed visually in Figs. P - §■ Plots showing 
the particle positions projected onto the equatorial plane are shown for Run 1 with n = 1.5 
in Fig. Run 2 with n = 1.0 in Fig. || and Run 3 with n = 0.5 in Fig. |7|. Corresponding 
contour plots covering the same spatial area in the equatorial plane are shown for Run 1 in 
Fig. |], Run 2 in Fig. [|, and Run 3 in Fig. || Time is measured in units of the dynamical time 
tj) for a spherical star with radius R eq as defined in Eq. ([11]); recall that these rapidly rotating 
models have significant rotational flattening. All models are rotating in the counterclockwise 
direction. Runs 1 and 2 were stopped at tf = 35£d and U = 50to, respectively, by which 
times these models had essentially stopped evolving. Run 3 was still evolving at if = 60£d 
when it was stopped because of the need to save computer time. Table [II] displays several 
important parameters of these models. 



The three models exhibit certain basic features in common; c.f. ||i3| , p^ , |2Q , |2^ . [2^ , |26| . Non- 
axisymmetric structure grows spontaneously out of deviations from axisymmetry, caused 
in these models by particle discreteness. A global bar-shaped structure develops in which 
the amplitude of the m = 2 mode grows exponentially in time; see Sec. [V]. During the 
bar mode's growth, trailing spiral arms develop as mass is shed from both ends of the bar. 
The bar and spiral arms exert gravitational torques that cause angular momentum to be 
transported outward from the core and lost from the ends of the bar. The spiral arms 
expand supersonically, causing shock heating and dissipation. Careful examination of the 
contour plots shows that the cores recontract toward an axisymmetric state after the initial 
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growth of the bar and ejection of spiral arms. These systems undergo two or more such 
episodes, depending on n; we will describe this in more detail below. Table [TTT] shows that 
the amounts of mass and angular momentum in the cores at the ends of the Runs 1 and 2 
are very similar, as are the final values of the stability parameter (3. Here we define the core 
to contain material within cylindrical radius w = R eq , where R eq is the initial equatorial 
radius. Since Run 3 was still evolving when the simulation was ended, we do not know what 
its final values will be. Overall, a significant amount of angular momentum is removed from 
the core by a relatively small amount of mass. 

Figure || displays the mass and angular momentum distributions for the three runs in the 
initial [frames (a) and (b)] and final [frames (c) and (d)] states. Here, m{w)/M and J{w)/ Jq 
are the normalized mass and angular momentum within cylindrical radius w, respectively. 
M is the total mass and Jo is the initial total angular momentum. Notice that at the 
final times the curves for both the mass and angular momentum distributions intersect near 
w ~ R cq . This is consistent with the fact that the models all shed roughly the same total 
amount of mass and angular momentum, and have final cores with radii w ~ R eq . 

Figs, ^(l) and ^(u) show that the final states of both Runs 1 and 2 exhibit a flattened 
"double halo" structure, consisting of a denser, inner region surrounded by a more diffuse, 
extended outer distribution of matter. (We do not know if Run 3 will have the same double 
halo structure at the conclusion of its evolution since the model was not run long enough to 
reach its final state.) This double halo may result from differences in the angular momentum 
carried by the mass when it is shed. During the first spiral arm ejection phase, the system 
has a higher value of (3 and the mass is shed from the ends of the bar with a greater angular 
momentum. This mass moves into the vacuum carrying a fraction of the initial angular 
momentum of the system, and eventually distributes itself about the remaining core to 
produce the outer halo. The inner halo is formed when the system undergoes the second 
spiral arm ejection episode at a lower value of /3, and the mass is shed with lower angular 
momentum into a region that already has mass from the first episode. 



For example, Fig. |10| shows the mass Am/M and angular momentum A J / Jq distributions 
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for Run 1 at t = 13.4£d (solid line) and t = 18.5^ (dashed line); these times correspond 
to the first and second episodes of mass shedding through spiral arms and are shown in the 
contour frames (e) and (i) in Fig. |j. Here Am is the amount of mass within a cylindrical 
shell of thickness dm at radius m, and similarly for A J. Fig. |i~0|(a) and (b) show the first set 
of spiral arms represented as a localized concentration of mass directly outside the stellar 
core in the region 1 < m/R cq < 2 at t = lSAto- As the model evolves, this mass expands 
into the surrounding vacuum. The amount of mass shed during the second episode is much 
less than during the first; c.f. Table [IV|. This can be easily seen when we examine Fig. |iO|(b), 



which zooms in on the ejected mass. Fig. |lO|(c) and (d) show that the material ejected from 
the primary instability carries a higher amount of angular momentum than the mass ejected 
after the second spiral arm ejection phase. 

Comparison of Figs. [3] - |8| shows that several important model properties depend on 
the stiffness of the equation of state. For example, the spatial deformation of the initially 
axisymmetric model into an elongated bar-shaped figure increases as n decreases and the 
fluid description approaches that of an incompressible fluid pi|. The widths of the spiral 



arms and bar also depend on the polytropic index, both decreasing as n decreases. And, as 
already mentioned, the system undergoes more spiral arm ejection phases as the equation 
of state stiffens. 

Recall that for stiffer polytropes, the density profiles become less centrally condensed, as 
shown in Fig. [I]. This greater amount of mass near the stellar boundary causes the material 
at the edge to be more tightly bound. Also, since the same rotation law Eq. is applied to 
all runs, the less compressible models exhibit a smaller degree of differential rotation. This 



can be seen by examining Figs. |2] and [11], which show, respectively, the angular velocities 
of the models at the initial and final times. The results of these effects can be clearly 
seen in frame (c) of Figs. |3], |5], and [7], which are the particle position plots for Runs 1, 2 
and 3, respectively. Comparison of these frames shows the mass at the stellar boundary 
becoming less diffuse as the polytropic index decreases. Overall, these effects contribute to 
the development of the less tightly wound spiral arm pattern found in models with lower 
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polytropic indices [pit . 



The behavior of the stability parameter for these runs is shown in Figure [12]. The solid 
lines give (3 = T TOt /\W\ and the dashed lines show T/|W|, where T rot is the rotational kinetic 
energy and T is the total kinetic energy Comparison with Figs. - [8] shows that (3 decreases 
sharply from its initial value as the bar instability develops, dropping to a local minimum 
as the bar reaches its maximum elongation. For stiffer polytropes, the temporal location 
of the initial growth of the bar occurs later in the evolution, with the minimum value of (3 
occurring at ~ 13to, 17to, and 19io for Runs 1, 2, and 3, respectively. When the results 
of Williams and Tohline |26| are converted to time measured in units of the dynamical time 



£dj they also show that the instability reaches a nonlinear amplitude later as n is decreased. 
Notice also that the value of the first local minimum of (3, which corresponds to the end of 



the first spiral arm ejection phase, decreases with n; c.f. |pq| . This behavior reflects the fact 
that as n decreases, the maximum elongation of the central bar-like region, and hence its 
moment of inertia, increases. Assuming angular momentum conservation in the core, this 
causes the minimum kinetic energy to decrease with n. 

At the end of the first spiral arm ejection phase, the core recontracts and (3 increases 
again. As the models evolve forward in time, they undergo successive periods of spiral arm 
ejection and core recontraction. The number of these episodes increases as the equation of 
state stiffens, with Run 1 showing 2 spiral arm ejection phases and Run 2 showing 4. Run 
3 undergoes 5 such episodes before the run was stopped; we expect it would exhibit several 
more if allowed to run to later times. Table [IV| shows the cumulative amount of mass and 
angular momentum shed from the core {w < R eq ) after each spiral arm ejection phase. 
Notice that the cumulative amount of angular momentum lost by the core after each spiral 
arm ejection episode decreases as the equation of state stiffens. Therefore, assuming that the 
cores conserve angular momentum between periods of spiral arm ejection, the stiffer cores 
recontract to a higher angular velocity (and a larger (3), and deform to a greater elongation 
(and a smaller (3) than more compressible ones, as shown in Fig. |T^. Also, since the cores 
lose angular momentum with each spiral arm ejection, each successive episode produces a 
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smaller maximum and a larger minimum value of (3. 

Now consider Run 3, which undergoes 5 periods of core recontraction corresponding 
to the local maxima of (3 seen in Fig. |12[ Its contracted core, displayed in the contour 
plots in frames (g), (k), (o), (r), and (w) of Fig. |8|, shows a "parallelogram-like" structure. 
This feature becomes stronger as n is decreased, with an intercomparison of frame (g) of 
Figs. f|, H and [| showing the emergence of this feature as we progress through the polytropic 
sequence toward stiffer equations of state. An explanation of this feature is as follows. When 
core recontraction begins, the ends of the bar move toward the central regions. A more 
compressible model is better able to increase its central density in response to this material 
forced toward the center. However, as n decreases, the material in the center cannot be easily 
compressed, thus forcing the fluid to move in a direction perpendicular to the contracting 
bar. This produces the observed parallelogram-like structure, or "anti-bar". We shall see in 
Sec. [V] below that the m = 4 mode is also present in these simulations; this may provide a 
degree of freedom that allows the formation of this feature. 

Another important difference observed as n is changed concerns the long term behavior 
of the models. The spiral arms in Runs 1 and 2 eventually merge as the systems evolve, 
resulting in a late time state consisting of a nearly axisymmetric central remnant of extent 
~ -R eq surrounded by a flattened double halo. After a comparable period of time, the core 
of Run 3 was still quite elongated and the spiral arms were just beginning to merge. One 
explanation for this longer-lived elongation in the n = 0.5 case is as follows. Consider an 
equilibrium sequence of uniformly rotating axisymmetric polytropes parametrized by (3. As 
(3 increases along such a sequence, a point is eventually reached at which mass is lost at the 
equator. Uniformly rotating polytropes with n > 0.808 (r < 2.24) reach this mass-shedding 
limit before the point at which ellipsoidal configurations can exist |pTT]JT8|] . Although Fig. |TT 



shows that the central remnants in these runs are differentially rotating, we believe that a 
similar mechanism may be operating here (see also p9|j4*0[] ), causing the cores of Runs 1 and 
2 to be nearly axisymmetric at the end of the run. 

One major difference observed between our work and that of Tohline, Durisen and collab- 
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orators is the final outcome of the simulations. In all of our runs with n = 1.5 and n = 1.0, 
we find that the systems evolve to a final state consisting of a nearly axisymmetric central 



remnant surrounded by an extended disk-like halo [2J),22]. This behavior was observed in 
both the n = 1.5 Eulerian and SPH runs investigated in the comparison study reported in 
Ref. [^3J as well as the work reported here. In contrast, all the long Eulerian evolutions 
reported by these other researchers (refs. |T7|Jl9|j20|j2^| ) resulted in a bar-like central core 
surrounded by a ring of material. (Interestingly, in the very low resolution SPH run re- 
ported in Ref. the low-density material did form an extended disk.) Currently, we do 
not understand the reason for these differences in the final state. One possible explanation 
is that these other researchers do not evolve an energy equation, and hence cannot model 



the shocks which occur in the outer, low density regions fl7| . Clearly, this is an important 
issue and efforts are underway to resolve these questions. 

V. ANALYSIS OF FOURIER COMPONENTS 

We quantify the development of the dynamical instability by studying the behavior 
of various Fourier components in the density using cylindrical coordinates (w,4>,z). The 
density in a ring of fixed w and z is analyzed using the complex Fourier series 

p{w,<j>,z) = C rn (zu,z)e m *. (15) 

m=— oo 

The azimuthal Fourier decomposition of the density distribution for various components m 
is expressed in terms of the amplitudes C m , defined by 

C m (w, z) = ±- r p(w, 0, z)e- m ^d<f). (16) 



2tt Jo 



4l| , |T9"|| . The relative normalized amplitude is then defined by 

\A m \ = \C m \/\C \, (17) 

where Co(m, z) = p(to, z) is the mean density in the ring under examination. The integration 
is performed over the azimuthal coordinate (0 < <p < 2tt) while w and z remain fixed. In 
this way, the analysis can be carried out in "density rings" for different values of w and z. 
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To apply this procedure to the SPH simulations, the particle model is interpolated onto a 
cylindrical grid at pre-chosen time intervals (typically every O.Olto) using kernel estimation 
p9| . The grid used here consists of 66, 34, and 16 zones in the zu, 4>, and z directions, 
respectively. Analysis of the density in rings at different values of vo in the equatorial plane 
(z = 0) gives quantitative information about the development of the global m = 2 mode 
visually seen in Figures § - ||, as well as other Fourier components that may be present. 

Examining how the normalized amplitude changes in time yields the growth rate 
d In | A m \/dt of the various Fourier components in our models. In practice, this is obtained by 
fitting a straight line through the data points in the time interval during which the function 
In \A m \ is linearly growing (thus giving an exponential growth rate for |A m |). The endpoints 
of this time interval are chosen "by eye" . A clearly defined linear region typically lasts for a 
relatively short time interval, and the value of the slope is sensitive to the endpoints defining 
the interval. 

Also, by examining the complex phase <j) m of a Fourier component, where 

Im(-C m 



(w, z) = tan 1 



(18) 



we can describe global non-axisymmetric structure propagating in the azimuthal direction. 
The development out of the initial noise of such a global mode with a well-defined angular 
eigenfrequency allows us to write 

4>m = cr m t, (19) 

where m is the phase angle of the disturbance, and a m is the eigenfrequency. The relation 
between the pattern speed fi pa t,m of the m th structure and the phase angle m is then plf 



o / \ 1 d(j) m <J m /on\ 

i2pat,m( W > Z ) = TT = • ( 20 ) 

m at m 

Notice that for the m = 2 bar mode, the eigenfrequency cr 2 is twice the rotational speed of 
the bar, and the bar rotation period is Tb ar = 2iim/a m = Att/^- 

The eigenfrequency is thus obtained by a simple calculation once the period is known. 
The period T m of the m th disturbance is determined directly from the period of the cosine 



of the phase angle (b m versus time. We use the function cos (b m rather than (b m itself due to 
the multi- valued nature of the inverse trigonometric function arctan (Eq. (|I8D), which would 
require us to artificially insert multiples of 7r in order to keep the function continuous. In 
practice, the half-period is obtained by locating successive differences in time between pairs 
of neighboring extrema of cos0 m . Once the half-period is known, the eigenfrequency a m can 
be obtained from Eq. (P0| ). Overall, we find that the m = 2 instability grows on a time scale 
of approximately one bar rotation period, as expected. 

The linearized tensor viral analysis (TVE) can also be used to calculate the bar mode am- 
plitude \A m \ and phase angle <p m . Although this method is exact only for small oscillations of 
uniform density, incompressible ellipsoids j|2]| , it has proven a useful point of comparison for 



numerical simulations when adapted to the study of rotating compressible fluids ||iT| , |H||2T 
Table [V] shows the TVE growth rates c? In | y4 2 1 / eft and eigenfrequencies 02 for the m = 2 
mode for the cases n = 1.5 and n = 1.0 with (3 = 0.31 reported by Williams and Tohline in 



Ref. PHI , where we have converted from their units. Notice that as n decreases, both the 
growth rate and eigenfrequency also decrease. We were unable to find TVE results in the 
literature for n = 0.5. 

Fig. [IB] shows the amplitudes of the Fourier components m = 1,2, 3, and 4 for Run 1 
in the density ring w = 0.36_R e q in the equatorial plane z = as functions of time. As 
expected, the development of the m = 2 disturbance dominates the initial evolution, with 
the other components growing at later times. Both the m = 2 and m = 4 components show 
an initial period of exponential growth. Since this takes place at various cylindrical radii w 
throughout the model, we identify these disturbances as global modes. The initial peak in 
the m = 2 amplitude corresponds to the maximum elongation of the bar and the minimum 
value of (3. The detailed structure after the initial growth of the bar mode varies somewhat 
with w, as the density in various parts of the star fluctuates due to the complex motions 
involved in the contraction and re-expansion of the core. 

The growth rate calculated for the m = 2 and m = 4 modes is rather sensitive to the 
specific time interval over which the linear fit to In \A m \ is performed. For the m = 2 mode in 
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Run 1, we typically get d\n \A^\/dt ~ 0.6^, and <iln l^l/dt ~ 0.8— li^ . The calculation of 
the eigenfrequencies is more robust, yielding cr 2 ~ lM^ 1 and <j 4 ~ S.Stf^ 1 . Both modes reach 
their peak amplitudes at about the same time, then drop to local minima and grow again. 
Since the pattern speeds of the these modes are nearly the same, f2 pati2 ~ fi pa t,4 ~ 0.95, this 
suggests that the m = 4 mode is a harmonic of the bar mode and not an independent mode 
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We also see that the m = 1 and m = 3 Fourier components grow somewhat, although 
not in the global and coherent fashion exhibited by the m = 2 and m = 4 modes. Recent 



work in the area of star formation [|T7i f43l has highlighted the importance of the ra — 1 case. 

Fig. [14| displays the amplitudes of the Fourier components m — 1,2, 3, and 4 for Run 
2 at the same value of w used above. Again, the m = 2 and m = 4 components emerge 
as exponentially growing global modes, with the bar mode dominating the early stages of 
the evolution. The growth rates in this case are even more sensitive to the time interval 
chosen for the linear fit than was the case for Run 1. We find <iln l^l/cfa ~ 0.5 — O.^^ 1 
and din \A 4 \/dt ~ 0.9 - 1.3^. A gain, the eigenfrequencies are less dependent on the time 
interval chosen and take the values o-i ~ l.St^ 1 and 04 ~ 2.9t5 1 - The pattern speeds are 
thus fipat,2 ~ ^ P at,4 ~ O.Tt^ 1 , implying again that the m = 4 mode is a harmonic of the bar 
mode. 

Finally, the amplitudes of the Fourier components m = 1, 2, 3, and 4 for Run 3 are shown 
in Fig. pl| The m — 1 and m = 3 Fourier components are stronger in this case, although 
they do not appear to develop into global modes. The m = 2 and m = 4 disturbances 
do develop into global modes and appear to be more strongly coupled than before. For 
example, the initial exponential growth rate of the bar mode (in the time interval 11.5^ < 
t ^ IStj} 1 ) is dh\\A2\/ dt ~ 0.4t5 1 - Then, the bar mode growth rate increases sharply to 
d\n \ A2\/dt ~ 2.2t^; this may due to coupling with the m = 4 mode, which initially grows 
at the rate din ~ 2.2t5 1 - The eigenfrequencies are 02 ~ 1.2^ and a 4 ~ 2.3^, so 

that fi pat ,2 ~ fi P at,4 ~ OM^ 1 - 
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Overall, the amount of structure seen in the Fourier components increases as n decreases. 
This reflects the fact that the stiffer fluids show more internal fluctuations as the cores expand 
and recontract. The eigenfrequencies a "2 do show the decrease with n predicted by the TVE 
analysis as given in Table |V|. This is due to the fact that the stiffer equations of state 
produce longer bars, which rotate more slowly. However, it is more difficult to assess the 
trends in the growth rates of the bar mode. If we consider the initial exponential growth 
period of the bar mode in Run 3, then it does grow at a slower rate than in Run 1 until 
the m = 4 mode starts to grow. As the fluids become stiffer and the number of spiral arm 
ejection episodes increases, the matter in the cores oscillates more. The coupling between 
the m = 2 and m = 4 also grows stronger; this may be linked to the development of the 



anti-bar discussed in Sec. IV. In particular, the elongation of the anti-bar in Run 3 seen in 



Fig. |8] (g) at t = 25.2£d occurs at roughly the same time as the second maximum in In \A±\ 



shown in Fig. 15 



VI. GRAVITATIONAL RADIATION 

The time-changing quadrupole moment caused by the development of the bar instability 
generates gravitational waves. The initial development of the bar mode produces a burst of 
radiation, followed by a weaker signal due to the subsequent expansions and recontractions 
of the core. Overall, the gravitational wave signal lasts for a longer time as the equation of 
state stiffens and the systems undergo more episodes of spiral arm ejection. Some interesting 
properties of the gravitational radiation produced by these models are given in Table [V]]. 

The gravitational waveform rh + for an observer on the axis at 9 = 0, = of a spherical 



coordinate system centered on the source is shown in Fig. [16] for these runs. Comparison 
of the waveforms with the contour plots in Figs. f|, || and [8] shows that indeed the onset 
of the burst coincides with the development of the primary instability. Notice that the 
maximum amplitude of the waveform does not vary significantly with the equation of state; 
see Table [V]]. As the core recontracts back to a more axisymmetric state, the amplitude of 
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the waveform decreases. The successive periods of spiral arm ejection and core recontraction 
produce additional bursts of gravitational waves; these have decreasing amplitudes because 
the maximum elongation of the core drops with each episode. Runs 1 and 2 show the weak 
signal of a slightly non-axisymmetric remnant in their final states, whereas Run 3 shows the 
much stronger signal of its still-evolving, more elongated core. Also, as the equation of state 
stiffens the frequency of the waves decreases, since the more elongated bars produced for 
smaller n rotate more slowly. 



The gravitational wave luminosity L for these runs is displayed in Fig. [jj]. Notice that 
the peak amplitude of the luminosity decreases as n decreases. For a non-axisymmetric 
object rotating rigidly about the z axis, the luminosity takes the form 

L = d 4 =-!?<'■-'.>"•■ (2D 

where I x and I y are the moments of inertia about the x and y axes, respectively, and Q is the 
rotational angular velocity ]I2 . Since the term in Q 6 dominates, the peak luminosity should 



decrease as n decreases and the central bar-like structures rotate more slowly, as shown in 
Fig.0 

It is interesting to examine the structure of the luminosity profiles. The luminosity 
of Run 1 shows peaks at t ~ 13£d an d t ~ 19£d! these correspond to the primary and 
secondary spiral arm ejection episodes. Run 2 shows two closely spaced peaks at t ~ 15to 
and t ~ 19tD, followed by other peaks at t ~ 25£d and t ~ 33£d- The first two peaks are 
associated with the initial period of spiral arm ejection; c.f. Fig. [12] (b). The remaining two 
peaks correspond to subsequent episodes. The successively smaller amplitudes reflect the 
fact that the angular velocity of the core decreases as angular momentum is shed on each 



subsequent episode. Finally, Fig. (c) shows that Run 3 has two closely spaced luminosity 
peaks at t ~ 17trj and t ~ 2lt-£>, which are again associated with the initial burst. The local 
miniumum in the luminosity at t ~ 25£d occurs at the time of core recontraction, as can be 
seen by comparing with Figs. |S| (g) and |T2| (c). At later times it is more difficult to discern 
individual bursts in the luminosity function, a trend that is also seen in the waveform shown 
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in Fig. H (c). 

The energy emitted as gravitational waves AE/Mc 2 is shown in Fig. [18]. In the case of 
Runs 1 and 2, AE/Mc 2 grows due to the initial and secondary bursts, and levels off when 
the cores reach their nearly axisymmetric final states. For Run 3, this quantity grows almost 
linearly with time and has not yet leveled off by the end of the simulation, indicating that 
the core is still quite nonaxisymmetric. 



Fig. [19] shows the rate at which angular momentum is carried by the waves, dJ z /dt. 
As was the case with the luminosity, we see structure in this quantity that corresponds to 
periods of spiral arm ejection and core recontraction. The angular momentum A J z carried 



by the gravitational waves is displayed in Fig. |2C] and shows features similar to those found 
in AE/M. 



Finally, the gravitational wave energy spectrum dE/df is displayed in Fig. ^Tjand shows 
that the peak frequency of the gravitational radiation / grav decreases as the equation of state 
stiffens. Table |V11| shows that, as expected, 2/ bar ~ / gr av, where /bar = (l/2)a 2 /27r is the 
rotational frequency of the bar. Notice, however, that the rotational frequencies 2 /bar are 
slightly lower than f gTllv . This is due to the fact that, while the the eigenfrequency cr 2 is 
calculated only during the initial development of the bar instability, / grav is computed for 
the entire evolution of the model and thus includes the higher rotational velocities obtained 
when the cores recontract. 



VII. SUMMARY AND DISCUSSION 

We have carried out numerical simulations of the dynamical instability in rapidly rotating 
stars initially modeled as polytropes with n = 1.5, 1.0, and 0.5. These calculations have 
been done using a 3-D SPH code with iV ~ 16, 000 particles. The code has a purely 
Newtonian gravitational field, and the gravitational radiation is calculated in the quadrupole 
approximation. The back reaction of the gravitational radiation is not included. 

All models exhibit the growth of the global m = 2 bar mode, with mass and angular 
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momentum being shed from the ends of the bar to form two trailing spiral arms. In general, 
as n decreases the central bar becomes narrower and more elongated. Once the central core 
has reached its maximum elongation, it begins to recontract toward a more axisymmetric 
state. This primary instability is followed by successive episodes of spiral arm ejection and 
core recontraction, with the number of these episodes increasing for stiffer equations of state. 
At the end of the simulations, the models with n = 1.5 and n = 1.0 have settled into a state 
with a nearly axisymmetric core of radius ~ R eq , where R eq is the initial equatorial radius, 
surrounded by a flattened disk-like halo that contains ~ 10% of the total mass and ~ 30% 
of the total angular momentum. Since these models have (3 S < (3 < (3d, they are expected 
to continue evolving under the secular instability pHf . The model with n = 0.5 had a fairly 



elongated core and was still evolving when that run was terminated. 

The development of the instability produces a burst of gravitational radiation. The 
maximum amplitude of the waveform r\h\ does not vary significantly with the polytropic 
index, whereas the frequency of the waves decreases somewhat as n decreases. This lowering 
of the frequency with n reflects the fact that the stiffer polytropes produce more elongated 
bars, which rotate more slowly; it also results in a decrease in the peak gravitational wave 
luminosity with n. Since the stiffer models undergo more episodes of spiral arm ejection and 
core recontraction, they produce longer-lived gravitational wave signals from the dynamical 
instability, with the total amount of energy and angular momentum emitted in the form of 
gravitational radiation increasing as n decreases. The nearly axisymmetric final cores (for 
n = 1.5 and n = 1.0) will continue to emit gravitational radiation as they evolve under the 
secular instability; this has been calculated by Lai and Shapiro [f|4] . 



The actual values of the gravitational wave quantities depend sensitively on the equatorial 
radius R eq of the stellar core when the dynamical instability takes place. This in turn depends 
on the astrophysical scenario in which the instability develops. Consider, for example, the 
collapse of a rotating stellar core of mass M = 1.4M that has (3 > (3d and is prevented 
from collapsing further due to centrifugal forces. The equatorial radius R eq of the core at 
which this centrifugal hangup occurs determines the amplitude and frequency of the resulting 
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gravitational radiation. The simulations presented here use stiff equations of state, which 
are appropriate only for stellar cores that have collapsed to near neutron star densities. We 
therefore calculated the gravitational wave amplitudes and frequencies from our models for 
two representative values of this parameter, R eq = 10 km and R eq = 20 km. We remind the 
reader that since these simulations have been done in the Newtonian limit, which breaks 
down for R eq ~ 10 — 20 km, these results must be viewed with appropriate caution. 

Table |VIII| shows the maximum amplitudes \h\ of the gravitational waveforms and the 



characteristic frequencies / grav for these representative values of R eq . Wave amplitudes are 
given for sources within the Milky Way (r = 15 kpc), the Local Group (r = 1 Mpc), and 
the Virgo Cluster (r = 20 Mpc). If the dynamical instability occurs at -R eq ~ 20 km, which 
is about twice the typical neutron star radius, / grav lies just outside the frequency range of 
the broad-band interferometers However, if such objects exist, they may potentially be 
observed using specially designed narrow-band interferometers |]45| , |46[| or resonant detectors 
jUJ . Of course, if hangup occurs at about the typical neutron star radius of R cq ~ 10 km, the 
characteristic frequencies become much larger. Since the star must be rotating differentially 
to achieve (3 > /3d, this last scenario could only occur in a newly- formed neutron star before 
its rotation becomes uniform (cf. ||12||). 

The maximum luminosity L/Lq, the energy emitted as gravitational radiation AE/Mc 2 , 
and the angular momentum carried by the waves A J / J is given for a core with mass 
M = 1.4M and these same representative values of R eq in Table |IX|. Here, L = c 5 /G 



and Jo is the initial total angular momentum. The largest integrated energy and angular 
momentum losses are produced by the model with n = 0.5. Since this model was still 
evolving when this run was stopped, the final values will be larger. 

There are several ways in which these calculations need to be improved to provide greater 
understanding of gravitational radiation from rotational instabilities. In this paper, we have 
concentrated on models with stiff equations of state. As noted above, these models are 
relevant for cores that have already collapsed to radii near the typical neutron star radius, 
R eq ~ 10 — 20 km. However, if centrifugal hangup occurs at R eq ~ 100 km the equation 
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of state is expected to be much softer, with n > 3. Simulations of this important case 
are currently in progress; these are being done using Eulerian techniques since we have 
found it easier to model the softer equations of state in this manner []4"7|j . Also, the very 
important and interesting question of the final state of the objects following the dynamical 
instability still remains to be fully resolved. In addition to longer runs with n = 0.5, this 
will involve a more detailed understanding of the differences between our results and those 
of Tohline, Durisen, and collaborators; we are making plans to pursue answers to these 
questions. Finally, gravitational radiation reaction and other general relativistic effects need 
to be included in order to have good physical models for comparison with future observations. 
We intend to include these effects in our future work. 



ACKNOWLEDGMENTS 

We are pleased to acknowledge interesting and helpful conversations with Steven Cran- 
mer, Dong Lai, Stephen McMillan, and Scott Smith. We are also grateful to Lars Hernquist 
for supplying a copy of TREESPH. This work was supported in part by NSF grant PHY- 
9208914. The numerical simulations were run at the Pittsburgh Supercomputing Center 
under grant PHY910018P. 



26 



REFERENCES 



[1] A. Abramovici, et al, Science 256, 325 (1992). 

[2] C. Bradaschia, et al, Nucl. Instrum. Methods A 289, 518 (1990). 

[3] K. Danzmann, et al, in Relativistic Gravity Research, Proceedings of the 81WE-Heraus- 
Seminar, Bad Hannef, Germany, edited by J. Ehlers and G. Schafer (Springer- Verlag, 
Berlin, 1992). 

[4] W. Johnson and S. Merkowitz, Phys. Rev. Lett. 70, 2367 (1993) 

[5] G. Harry, T. Stevenson, and H. Paik, Phys. Rev. D 54, 2409 (1996). 

[6] B. Schutz, in Dynamical Spacetimes and Numerical Relativity, edited by J. Centrella 
(Cambridge University Press, New York, 1986). 

[7] K. Thorne, in Compact Stars in Binaries, Proceedings of IAU Symposium 165, edited by 
J. van Paradijs, E. van den Heuvel, and E. Kuulkers (Kluwer, Dordrecht, 1996). 

[8] K. Thorne, in Proceedings of the Snowmass 95 Summer Study on Particle and Nuclear 
Astrophysics and Cosmology, eds. E. W. Kolb and R. Peccei (World Scientific, Singa- 
pore, 1995); also published in Particle Physics, Astrophysics & Cosmology, Proceedings 
of the SLAC Summer Institute on Particle Physics, eds. Jennifer Chan & Lilian DePorcel 
(SLAC-Report-484, Stanford Linear Accelerator Center, Stanford, CA, 1996). 

[9] R. Wagoner, Astrophys. J. 278, 345 (1984). 

[10] B. Schutz, Class. Quantum Gravity 6, 1761 (1989). 

[11] J. Tassoul, Theory of Rotating Stars (Princeton University Press, Princeton, 1978). 

[12] S. Shapiro and S. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars (Wiley, New- 
York, 1983). 

[13] R. Durisen and J. Tohline, in Protostars and Planets II, edited by D. Black and M. 



27 



Matthews (Univ. of Arizona Press, Tucson, 1985). 
[14] R. Managan, R. Astrophys. J. 294, 463 (1985). 

[15] J. Imamura, J. Friedman, and R. Durisen, Astrophys. J. 294, 474 (1985). 

[16] J. Imamura, J. Toman, R. Durisen, B. Pickett, and S. Yang, Astrophys. J. 444, 363 
(1995). 

[17] B. Pickett, R. Durisen, and G. Davis, Astrophys. J. 458, 714 (1996). 

[18] D. Lai, F. Rasio, and S. Shapiro, Astrophys. J. Suppl. 88, 205 (1993). 

[19] J. Tohline, R. Durisen, and M. McCollough, Astrophys. J. 298, 234 (1985). 

[20] R. Durisen, R. Gingold, J. Tohline, and A. Boss, Astrophys. J. 305, 281 (1986). 

[21] H. Williams and J. Tohline, Astrophys. J. 315, 594 (1987). 

[22] J. Houser, J. Centrella, and S. Smith, Phys. Rev. Lett. 72, 1314 (1994). 

[23] K. New, Ph.D. Thesis, Louisiana State University (1996). 

[24] S. Smith, J. Houser, and J. Centrella, Astrophys. J. 458, 236 (1996). 

[25] J. Houser, Ph.D. Thesis, Drexel University (1996). 

[26] H. Williams and J. Tohline, Astrophys. J. 334, 449 (1988). 

[27] J. Centrella and S. McMillan, Astrophys. J. 416, 719 (1993). 

[28] L. Lucy, Astron. J. 82, 1013 (1977); R. Gingold and J. Monaghan, Mon. Not. R. Astron. 
Soc. 181, 375 (1977); see J. Monaghan, Ann. Rev. Astron. Astrophys. 30, 543 (1992) for 
a review. 

[29] L. Hernquist and N. Katz, Astrophys. J. Suppl. 70, 419 (1989). 

[30] J. Barnes and P. Hut, Nature 324, 446 (1986); L. Hernquist, Astrophys. J. Suppl. 64, 
715 (1987). 

28 



[31] C. Misner, K. Thorne, and J. Wheeler Gravitation (Freeman, New York, 1973). 

[32] K. Thorne, in 300 Years of Gravitation, edited by S. Hawking and W. Israel (Cambridge 
University Press, New York, 1987). 

[33] L. S. Finn and C. Evans, Astrophys. J. 351, 588 (1990). 

[34] S. Smith and J. Centrella, in Approaches to Numerical Relativity, edited by R. d'Inverno 
(Cambridge University Press, New York, 1992). 

[35] J. Ostriker and J. Mark, Astrophys. J. 151, 1075 (1968). 

[36] P. Bodenheimer and J. Ostriker, Astrophys. J. 180, 159 (1973). 

[37] I. Hachisu, Astrophys. J. Suppl. 61, 479 (1986). 

[38] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes 
(Cambridge University Press, New York, 1992) 

[39] X. Zhuge, J. Centrella, and S. McMillan, Phys. Rev. D 50, 6247 (1994). 

[40] X. Zhuge, J. Centrella, and S. McMillan, Phys. Rev. D, submitted. 

[41] D. L. Powers, Boundary Value Problems, 3rd Edition, (Harcourt, Brace, and Jovanovich, 
Florida, 1987). 

[42] S. Chandrasekhar, Ellipsoidal Figures of Equilibrium, (Yale Univ. Press, New Haven, 
1969). 

[43] I. Bonnell, Mon. Not. R. Astr. Soc. 269, 837 (1994); I. Bonnell and M. Bate, Mon. Not. 
R. Astr. Soc. 271, 999 (1994). 

[44] D. Lai and S. Shapiro, Astrophys. J. 442, 259 (1995). 

[45] B. Meers, Phys. Rev. D 38, 2317 (1988); K. Strain and B. Meers, Phys. Rev. Lett. 66, 
1391 (1991); A. Krolak, J. Lobo, and B. Meers, Phys. Rev. D 43, 2470 (1991); 47, 2184 
(1993). 

29 



[46] D. Kennefick, D. Laurence, and K. Thorne, Phys. Rev. D, in preparation. 
[47] J. Centrella, J. Gao, and S. Smith, work in progress. 



30 



TABLES 



Model 


n 


N m 


N z 


Nt 


P 


R P 


/R eq 




VR 


SCF1 


1.5 


225 


225 


20 


0.30 


0.20 




7.6 x 10~ 5 


SCF2 


1.0 


129 


70 


16 


0.30 


0.23 




4.5 x 10~ 5 


SCF3 


0.5 


129 


70 


22 


0.30 


0.25 




4.0 x 10~ 5 


TABLE I. Properties of the initial axisymmetric equilibrium models created using the SCF 


method. 


N m and N 2 


are the number of uniform 


grid zones in the w 


and z 


directions, respectively. 


Nu is the number of iterations required for convergence to 


a solution with a 


tolerance of 10~ 5 . The 


axis ratio is R p /R eq 


The value of the virial parameter calculated on the c 


ylindrical grid is VR. 


Model 


n 


N 


Vi2|i 


A 


time 


E[ — i?f 




Ji-J { 


CPU 












[to] 








(hr) 


Run 1 


1.5 


16096 


4.5 x 10~ 2 


0.32 


35 


0.020 




< .001 18.8 


Run 2 


1.0 


16619 


3.7 x 10~ 2 


0.32 


50 


0.015 




.002 


25.2 


Run 3 


0.5 


16526 


1.5 x 10~ 4 


0.31 


60 


0.018 




.003 


29.0 



TABLE II. Properties of the SPH models. N is the total number of particles in each model. 
The fact that the method used to generate the initial particle models does not allow strict control 
over the the number of particles accepted into each model results in somewhat unusual values of 
N. The subscripts "i" and "f" denote the initial and final states of the model, respectively. The 
stability parameter of the particle model at the initial time is (3{. The duration of the run in units 
of the dynamical time to is given in the column labled "time" . E is the total energy, and J is the 
total angular momentum. All models were run on a Cray C90; the amount of CPU time used is 
given for the duration of the run. 
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Model 


±y - L core,! 


J r 

u core,! 


B f 

H core, r 


HI 




[%] 


[%] 






Run 1 


90 


70 


0.24 


0.26 


Run 2 


91 


70 


0.24 


0.25 


Run 3 


92 


71 


0.23 


0.24 



TABLE III. Hydro dynamical results for the models. The core refers to matter within cylindrical 
radius w = R cq , where R eq is the initial equatorial radius, and the subscript "f" denotes the final 
state of the model. 
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Run 


/ 

[to] 


Milled 
[%] 


<^shed 
[%] 


P 


Run 1 





0.0 


0.0 


0.32 




17 


7.3 


24 


0.27 




23 


8.6 


28 


0.26 




35 


9.5 


30 


0.26 


Run 2 





0.0 


0.0 


0.32 




22 


5.0 


18 


0.28 




30 


7.0 


25 


0.26 




36 


7.9 


27 


0.26 




42 


8.5 


29 


0.25 




50 


9.1 


30 


0.25 


Run 3 





0.0 


0.0 


0.31 




25 


4.1 


12 


0.29 




35 


6.2 


19 


0.28 




43 


6.7 


22 


0.26 




50 


7.1 


24 


0.26 




58 


7.4 


26 


0.25 




60 


8.4 


29 


0.24 



TABLE IV. Properties of the models after each successive spiral arm ejection phase. The mass 
-^shed and angular momentum J s hed shown are the cumulative mass and angular momentum lost 
after each such episode. The core is defined as mass within w = R cq ; see Fig. ||. The values of 
(3 are obtained directly from the successive peaks corresponding to core recontraction in Fig. [l^. 
The last temporal point in each series corresponds to the end of the run, and is not necessarily a 
time at which the core has reached maximum recontraction. 
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n 


din \A 2 \/dt tve 


0"2 TVE 








1.5 


0.73 


1.7 


1.0 


0.56 


1.5 



TABLE V. TVE bar mode growth rates d1n\A2\/dt and eigenfrequencies oi for n = 1.5 and 



n = 1.0 These values are taken from Ref. |21], where we have converted from their units 



Model max \rh\ max L/L (AE/Mc 2 )f max dJ z /dt (AJ z /J )f 

Run 1 0.57 0.13 0.87 0.066 1.02 

Run 2 0.58 0.091 1.1 0.057 1.53 

Run 3 0.58 0.059 2.2 0.044 3.41 



TABLE VI. Gravitational wave results for Runs 1, 2, and 3. The peak values of \rh\, L/Lq, and 
dJ z /dt throughout the run, and the final (cumulative) values of (AE/Mc 2 ) and (AJ Z / Jq) are given. 
Lq = c 5 /G and Jo is the initial total angular momentum. To obtain dimensional quantities, the 



scalings given in the axis labels of the corresponding Figs. [16| - |20| must be applied; see Tables [VTJ I 
and El. 
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Model 2/ bar / grav 

[tv] fc?] 

Run 1 0.30 0.32 

Run 2 0.24 0.30 

Run 3 0.19 0.23 



TABLE VII. Frequencies for the models, /bar is the rotational frequency of the bar and is cal- 
culated from the eigenfrequency oi- /grav is obtained from the gravitational wave energy spectrum 



dE/df shown in Fig. 21 





max \h\ww 


max \h lg 


max |/i|vc 


/grav 


/grav 


/grav 




(r = 15 kpc) 


(r = 1 Mpc) 


(r = 20 Mpc) 


(n = 1.5) 


(n = 1.0) 


(n = 0.5) 


10 km 


5 x 10~ 19 


8 x 10" 21 


4 x 10~ 22 


4900 Hz 


4100 Hz 


3100 Hz 


20 km 


3 x 10~ 19 


4 x 10~ 21 


2 x 10~ 22 


1700 Hz 


1400 Hz 


1100 Hz 



TABLE VIII. The maximum amplitudes of the gravitational waveform \h\ and the characteris- 
tic frequencies / gra v are given for two representative values of the equatorial radius R eq . The core is 
taken to have mass M = 1.4M Q . The waveform amplitudes \h\ are given for sources located within 
the Milky Way (r = 15 kpc), the Local Group (r = 1 Mpc), and the Virgo Cluster (r = 20 Mpc). 
Notice that \h\ is essentially independent of the polytropic index n. These values were obtained by 



applying the appropriate scalings to the data given in Tables VI and VII. 



35 





max L/Lq 


AE/Mc 2 


AJ/Jo 


10 km 


2 - 5 x 10~ 5 


4 - 9 x 10~ 3 


2 - 7 x 10~ 2 


20 km 


6 - 20 x 10" 7 


4 - 8 x 10~ 4 


4 - 10 x 10~ 3 



TABLE IX. The maximum luminosity L/Lq, the energy emitted as gravitational radiation 
AE/Mc 2 , and the angular momentum carried by the waves AJ/Jq are given for two representative 
values of the equatorial radius R cq . Lq = c 5 /G and Jo is the initial total angular momentum. The 
core is taken to have mass M = 1.4M Q . The lower and upper limits for L/Lq are produced for 
n = 0.5 and n = 1.5, respectively; the values for n = 1.0 are between these two limits. However, 
the lower values of AE/Mc 2 and AJ/Jq correspond to the case n = 1.5. The larger values are 
produced by the model with n = 0.5; since this model was still evolving when it was stopped, these 
values will be larger once the model reaches its final state. These values were obtained by applying 



the appropriate scalings to the data given in Table VI. 
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FIGURES 

FIG. 1. The normalized equatorial plane density is shown for the iterated contraction initial 
models Run 1 (n = 1.5), Run 2 (n = 1.0), and Run 3 (n = 0.5). Here, the equatorial plane is 
taken to include all particles within z = ±0.01R eq . The solid curve in each frame represents the 
SCF equatorial plane density, and p c is the the central SCF density. 

FIG. 2. The normalized angular velocity is shown for the initial models of Runs 1, 2 and 3. All 
particles are plotted in this figure. In each frame, the solid curve gives the angular velocity for the 
corresponding SCF initial model and Q c is the SCF central angular velocity. 

FIG. 3. Particle positions are shown projected onto the equatorial plane for various times 
during the evolution of Run 1 with n = 1.5. All particles are plotted. The vertical axis is y/R C q 
and the horizontal axis is x/R eq . The system rotates in the counterclockwise direction. 

FIG. 4. Density contours in the equatorial plane are shown for Run 1 with n = 1.5. The 
frames are taken at the same times as the corresponding particle plots in Fig. ||. The contour levels 
are the same in all frames, and are spaced a factor of 10 apart, going down 4 decades below the 
maximum (central) SCF initial density. The contours were calculated using kernel interpolation 
on a 100 x 100 Cartesian grid covering the same spatial area in the equatorial plane as the frames 
in Figure ||[ 

FIG. 5. Same as Fig. |3| for Run 2 with n = 1.0. 

FIG. 6. Same as Fig. ^ for Run 2 with n = 1. These frames correspond to the particle plots 
shown in Fig. |^. 

FIG. 7. Same as Fig. ||| for Run 3 with n = 0.5. 

FIG. 8. Same as Fig. ||| for Run 3 with n = 0.5. These frames correspond to the particle plots 
shown in Fig. |?j. 
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FIG. 9. The distributions of mass m{m)/M and angular momentum J{w)/Jq are shown for 
Run 1 (solid line), Run 2 (dashed line), and Run 3 (dot-dashed line). Frames (a) and (b) show the 
initial models, and frames (c) and (d) show the final states. Here, M is the total mass and Jo is 
the initial total angular momentum. 

FIG. 10. The mass Am/M and angular momentum AJ/Jq distributions are shown for Run 1. 
Here, Am is the mass within a cylindrical shell of thickness dw at radius w, and similarly for AJ; 
Jo is the total initial angular momentum. The solid lines show the values at time t = 13.4£d and 
the dashed lines at t = 18.5£d- Frames (b) and (d) show enlargements of (a) and (c), respectively. 

FIG. 11. The normalized angular velocity is shown for the final states of Runs 1, 2, and 3. 
Here, O c is the central angular velocity for the initial SCF model. All particles are plotted. 

FIG. 12. The behavior of the stability parameter j3 = T rot /| W| (solid line) and T tot /|VF| (dashed 
line) is shown as a function of time for Runs 1, 2, and 3. Here, T rot is the rotational kinetic energy, 
T to t is the total kinetic energy, and W is the gravitational potential energy. 

FIG. 13. The growth of the Fourier components m = 1,2,3, and 4 for Run 1 with n = 1.5. 
These values were obtained in the density ring at w = 0.36i? e q i n the equatorial plane z = 0. 



FIG. 14. Same as Fig. 13 for Run 2 with n = 1.0. 



FIG. 15. Same as Fig. |13| for Run 3 with n = 0.5. 



FIG. 16. The gravitational waveform r/i+ for an observer located at = <f> = at distance r 
from the source for Runs 1, 2, and 3. 

FIG. 17. The gravitational wave luminosity L/Lq for Runs 1, 2, and 3. Here, Lq = c 5 /G. 

FIG. 18. The energy AE/Mc 2 emitted as gravitational radiation for Runs 1, 2, and 3. 
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FIG. 19. The rate dJ z /dt at which angular momentum is carried away by gravitational radiation 
for Runs 1, 2, and 3. 

FIG. 20. The angular momentum AJ/Jq carried by the gravitational waves for Runs 1, 2, and 
3. Here, Jo is the total initial angular momentum. 

FIG. 21. The gravitational wave energy spectrum dE/df is shown as a function of frequency / 
for Runs 1, 2, and 3. 
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